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Abstract 

The process of condensation in the system of scalar Bosons with weak A</> 4 
interaction is considered. Boltzmann kinetic equation is solved numerically. 
Bose condensation proceeds in two stages: At the first stage condensate is 
still absent but there is non-zero inflow of particles towards p = and the 
distribution function at p = grows from finite values to infinity. At the 
second stage there are two components, condensate and particles, reaching 
their equilibrium values. We show that the evolution in both stages proceeds 
in a self - similar way and find the time needed for condensation, which is 
finite. 
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It is a fundamental result of quantum statistics of Bosons that above a certain critical 
density all added particles must enter the ground state: Bose-Einstein condensate forms. 
The kinetics of this process is a very interesting problem. One can reach the Bose condensa- 
tion gradually decreasing temperature in a sequence of equilibrium states. An appropriate 
description of this is given by the well-known kinetics of second order phase transitions. 
On the other hand, when the conditions for the formation of Bose condensate appear, the 
system can be far from the equilibrium. Recently this problem attracted particular interest 
in connection with exciting prospects for the experimental observation of Bose condensation 
in a very cold atomic samples, e.g. in a gas of spin-polarized atomic hydrogen |l] or in 
alkali- metal vapors @. Another interesting application of Bose kinetics is rather far from 
the laboratory experiments and is related to the problem of Bose stars formation from 
the dark matter in the universe. 

The question of the time evolution of weakly interacting Bose gas was addressed in a 
number of papers. In earlier treatments an ideal Bose gas was coupled to a thermal bath 
with infinite heat capacity |5||]]. Small energy exchange was assumed and after several 
other approximations a Fokker-Planck type equation was obtained. Levich and Yahkot 
H calculated analytically that the time for condensation is infinite in this situation. By 
including Boson-Boson interactions they later Q found solution which describes explosive 
appearance of a condensate, but they concluded that this effect could have been an artifact 
of their approximations. 

Snoke and Wolfe in Ref. || undertook direct numerical integration of Boltzmann ki- 
netic equation. Although this calculation demonstrated the restructuring of the distribution 
function, the appearance of a Bose condensate was not detected. Their approach is close to 
ours but in comparison with the Ref. || we perform numerical integration in much wider 
dynamical range of relative energies and densities and we directly analyse the behavior the 
distribution function. 

Another result recently reported in Ref. states that the time required for condensation 
is ~ T c , where T c is the temperature of Bose condensation. This result seems to be incorrect 
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(or at least it can not be applied to all temporal stages of gas evolution) since it is insensitive 
to the interactions. Clearly, in the limit of zero couplings the distribution of particles does 
not evolve and the relaxation time has to be infinite. 



Analytical study of Bose condensation was performed recently in the paper [K| . Three 
different regimes of evolution where identified. It was argued that in the kinetic region in 
non-linear regime the distribution function has to be a power law / oc e- 7/6 . This power low 



is well known in the theory of plasma turbulence [fn] , |12|| . We indeed observe the tendency 
to this law in our numerical simulations, but the system never reaches it. Moreover, this 
distribution became destroyed with condensate appearance, contrary to the assumptions of 
Ref. [[HJ. What concerns condensation time, only rough dimensional estimations were done 
in Ref. [K|. 

We found that Bose condensation process can be divided on two stages. During the 
first part of the first stage the self-similar solution forms and then the distribution function 
reaches infinity at p = in a self-similar way. We can not observe in principle the actual build 
up of the coherence in the frameworks of Boltzmann equation. But we have to conclude 
that after the distribution function became infinite at zero momentum, the condensate 
had formed. While due to fluctuations the width of the coherent region can be finite in 
momentum space, we simply model the coherent field by N c S(p) with initially infinitesimally 
small but then growing amplitude N c (t). We found that in this second stage the evolution 
also proceeds in a self-similar fashion and found the durations of both stages. 

While our prime interest and motivations for this work was connected to the physics of 
Bose star, the problem of condensate formation with account for the gravitational field was 
never approached and we do not attempt it here. Instead, we are solving kinetic equations 
in flat space-time and only range of parameters of the initial distribution reflect the virial 
equilibrium of self-gravitating system. 

We consider the system of scalar bosons with 4-particle self-interaction. The field- 
theoretical Lagrangian which we bear in mind is L = (<9 M 0) 2 /2 — m 2 2 /2 — A0 4 /4!. Time 
development of the quantum state which contain well defined and large number of particles 



can be adequately described by the Boltzmann kinetic equation which governs the evolution 
of one-particle distribution function, /(p): 



where 

F(/) = [1 + /x] [1 + f 2 ]f[f> - [1 + /{] [1 + /^]A/ 2 , (2) 

and fi = /($), = f{pi')- The equilibrium solution of the kinetic equation is the Bose- 
Einstein distribution function 

/ (P) = Tt \m i + ( 27r ) 3iV ^(P) > (3) 

exp ((e - n)/T) - 1 

where e is the particle energy, \i is the chemical potential, T is the temperature of the final 
state and N c is the number density of particles in condensate. 

In what follows we shall consider isotropic initial distribution f = f(e). In our case 
the matrix element is given by |M/j| = A 2 and the kinetic equation for the case without 
condensate, N c = 0, can be rewritten in the form 

where D = mm\pi, p 2 , p[, j/ 2 ] and e 2 = e[ + e' 2 — E\ in arguments of F(f). The integration 
should be done over the region e\ < e[ < oo, e x — e[ < e' 2 < oo. 

After the moment of condensate formation the kinetic equation (|]) is inappropriate for 
numerical integration anymore, the finite number of particles in condensate corresponds 
to the infinite value of the distribution function at zero energy. In order to describe the 
system of particles interacting with the condensate we divide the distribution function into 
two pieces: / = f(s,t) + (27r) 3 iV c (t)<5 3 (p), where the first term corresponds to the "gas" of 
particles and the second one describe the condensate. Substituting this function into the 
original kinetic equation ([!]) we obtain 
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N c (t) = r de' x de' 2 [f x f 2 - f 2 (f[ + %)} , (5a) 

bArr^m Jo 

f(ei) =I P + ( - h)A - fihW, + 2 - k)h + mde' 2 ) . (5b) 

In general, after the condensate formation (and at large particle densities even before) 
the description in terms of quasiparticles rather then particles is more appropriate. For 
example, the kinetic equation, Eq.@, does not include the processes where one of the 
incoming and one of the outcoming particles has zero momentum, p' 2 = p 2 = 0. This process 
does not contribute to the collision integral directly, i.e. it does not change the distribution 
of particles over energies, but it does change the effective particle mass. 

However, in many cases those effects are insignificant and we still can work in terms of 
particles. The quantitative arguments are the following. The effective mass of quasiparticles 
in the presence of condensate (or in dense medium) is m e fr = m 2 + \<p 2 /G = m 2 + Xn/m. 
We can still use the description in terms of particles if the second term in the sum is much 
smaller than the first one. Since n ~ m 3 (A-u) 3 /o we obtain A/o(^^) 3 ^ 1> where At; is 
characteristic velocity dispersion. In the case of axion miniclusters, for example, we have Q 
A/o ~ 10 -5 ; A-u ~ 10~ 8 , and the description in terms of particles is perfectly good. 

As an initial distribution we choose the function f(e) which has the maximum at p = 
0. In general, such distribution function can be characterized by means of three major 
parameters: (1) The overall amplitude /q. In what follows we define /o = f(e = 0). (2) The 
energy scale eo where the distribution function became twice smaller, f(eo) — fa/2. (3) The 
effective width, T, of the region over which the distribution function varies rapidly. More 
specifically, we choose the initial distribution function to be of the form: 

2 y 

f(e) = arctan [exp (r(l — e/s ))] . (6) 

7T 

In what follows we shall measure the energy in units of Eo and the distribution function in 
units of fo, i.e. initial distribution function will have the normalization f(e = 0) = 1. 
We define the dimensionless time r as 0: 
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e 2 foil + /o)A 2 



t . 



(7) 



647r 3 m 



The parameters Eo and fo in the limit /o ^> 1 after rescaling will not enter the kinetic 
equation explicitly, but will define the time scale. In terms of r there remains the weak 
dependence of the relaxation time upon the initial shape parameter T. All data presented 
in this paper will correspond to one and the same value of T = 5. With the initial shape 
of the distribution function being given, the two scaling parameters fo and £o define also 
the parameters of final equilibrium {T, /i, N c }. With T = 5 the choice fo < f cr it ~ 2.8 
corresponds to the Bose-gas without condensate, while fo > / cr it corresponds to condensate 
formation in the final equilibrium state. 

We shall consider here the case fo > 2.8 only. We shall simplify the problem and 
consider fo ^> 1. In this case we can disregard f 2 terms in the function (0), which became 



We integrated the kinetic equation in the energy interval 10~ 9 < e < 10. We defined the 
distribution function on the grid of 200 points equally spaced in the logarithm of energy and 
used the spline interpolation when calculating distribution function at intermediate points. 
For each integration in collision integral we had implied Gauss algorithm. Particle and 
energy non-conservation was of order 10~ 3 for the whole time of integration. 

Results of numerical integration of the kinetic equation @ are presented in Fig. 1, 
where we plot the distribution function at different moments of time. We have arranged 
the output each time when f(e = £ m i n ,t) had increased by one order of magnitude. The 
most striking feature of this plot is self-similar character of the evolution. The distribution 
function has the "core" where f(e) ~ const and the radius of "core" decreases with time 
while the value of /(e) in the core grows. Outside the core the distribution function is the 
power law /(e) oc E~ a and does not depend upon time to a very good accuracy. Self-similar 
solutions exhibiting this kind of behavior can be parametrized as 



p(f) = [fi + wvi - \f[ + mih. 



f(E,r)=A- a (r)f s (E/A(r)) , 



(8) 



where it is assumed A{r) 



with the increase of time and we always can choose the 
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normalization / s (0) = 1. It is possible to find the time dependence of distribution function 
at e = at late times analytically, using self-similarity of the solution. 
Substituting parametrization (H) in the kinetic equation we obtain: 

/(0,r) = [2C(r c -r)(«-l)]-^- 1 )) , (9) 

and f(0, t) has the pole at the finite moment of time r = r c , f(0, r c ) = oo. 

The value of logarithmic derivative of f(e) is plotted in Fig. lb. The dotted line in Figs. 
1 corresponds to the limiting value a = 7/6. This power law corresponds to the stationary 
solution of the kinetic equation (f|), see Ref. ||12|| . But this value is never reached prior to the 
moment of condensate formation (after that moment the character of evolution completely 
changes). Rather, with the boundary condition df/de = at e = the power law on the tail 
is a m 1.24 (with appropriate boundary conditions our code correctly finds the root a = 7/6 
). And the function Eq. @ is a good fit to our data with a ~ 1.24 and r c ~ 19. 

We conclude, that after the solution have reached self-similar form, the time dependence 
of the distribution function at zero momentum is given by /(0, r) oc (r c — t)~ 2 ' 6 , which 
reaches infinity at finite time and means the onset of condensation. 

Though the system does not reach the limiting value of a = 7/6, at least this power law 
is the root of the equation Ip = [12|]. Because of that we selected the function f(e) oc £~ 7 / 6 



at e < 1 as an initial condition while integrating Eqs. @. We took initially N c <C iV tot (the 
particular choice for N c , as far as this condition is satisfied, is insignificant). 

Results of numerical integration of the system @ are presented in Fig. 2. The dashed 
line corresponds to the initial distribution. Solid lines correspond to the distribution function 
at different moments of time. Basically, evolution proceeds in the following way. First, the 
power law f(e) oc l/(e) 7 ^ changes to the law f(e) oc \ je at small energies. And then this 
change propagates to the region of lager energies, see Fig. 2. Later on the power law stays 
on the equilibrium value a — 1, but the amplitude of f(e) gradually decreases. 

Again, we see that the essential parts of curves in Fig. 2 repeat itself under translation 
from the left to the right and evolution is self-similar. During this epoch (before "the wave of 



change" had reached exponential tail of the initial distribution at e > 1) approximately 40% 
of particles had condensed. And what is important, the number of particles in condensate 
linearly grows with time at this epoch, N c /N tot = Br. This enables us to eliminate the 
ambiguity in the initial value for N c since B does not depend upon it. Indeed, in our 
simulations which was done in a finite energy interval, during several first iterations system 
self- adjusts: a proper profile of distribution forms while condensate reaches particular value 
of N c . We can disregard this period and extrapolate curves in Fig. 2, N c (t) and self-similar 
character of the evolution back in time and to the region of smaller energies. 

We obtained B « 2 x fi(l) 2 (initial normalization fi(l) = 1 corresponds roughly to 
the magnitude of final distribution at e ~ e min in Fig. 1). For completeness we also had 
runs when we used final distribution in Fig. 1 as initial condition for integration Eqs. (|5|). 
Condensation in this case decelerates and N c does not grow lineraly, all is in accordance 
with excess magnitude of the distribution function at small energies when a ~ 1.24. 

We have studied numerically the kinetics of Bose condensation of the weakly interacting 
Bose gas. The picture we observe differs even qualitatively from the previous works. The 
distribution function of excess particles which eventually has to form the condensate does 
not narrows in time gradually approaching 5- function, as it was found in Ref. |J assuming 
small energy exchange per collision. Instead, power low profile / ~ e~ 7 ^ tend to form which 



corresponds to the constant flux of particles in momentum space towards the condensate [12 



One could expect, as it was done in Ref. ||10|| , that this very natural regime will persist in the 



presence of condensate as well, till all excess particles will inflow from the high energy tail 
and reside in the ground state. Nevertheless, this is not the case also and at the moment the 
condensate appears, it terminates this regime. Instead of this steady flow through the whole 
energy interval, particles from all energy levels jump directly to the condensate retaining 
during the major period of time the equilibrium shape of the distribution function, / ~ e^ 1 . 
The constant of proportionality in this law gradually decreases till it reaches equilibrium 
value. 

We can not, in principle, observe the build up of coherence being in frameworks of 
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Boltzmann kinetic equation, but kinetic description has to be valid prior to and after the 
very moment of condensate formation. We have seen that evolution on both stages is self- 
similar. This allows us to obtain a number of useful analytical relations, e.g. the time 
dependence of the distribution function near the point of condensate formation, Eq. []. We 
have found the duration of both stages, which is finite and of order r c ~ 20. 

We thank E. Kolb, V. Rubakov and G. Starkman for useful discussions. This work was 
supported by DOE and NASA grant NAGW-2381 at Fermilab. D.S. thanks the Astrophysics 
Department of FNAL for the hospitality during this work. 
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FIG. 1. The distribution function (top) and its logarithmic derivative (bottom) are shown at 
different moments of time prior to condensate formation. 
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FIG. 2. The distribution function (top) and its logarithmic derivative (bottom) are shown at 
different moments of time during condensation 
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